#define Ntuple_cxx
// The class definition in Ntuple.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.

// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// Root > T->Process("Ntuple.C")
// Root > T->Process("Ntuple.C","some options")
// Root > T->Process("Ntuple.C+")
//

#include "Ntuple.h"
#include <TH2D.h>
#include <TH1D.h>
#include <TStyle.h>

 void BinLogX(TH1*h)
{
   TAxis *axis = h->GetXaxis();
   int bins = axis->GetNbins();

   Axis_t from = axis->GetXmin();
   Axis_t to = axis->GetXmax();
   Axis_t width = (to - from) / bins;
   Axis_t *new_bins = new Axis_t[bins + 1];

   for (int i = 0; i <= bins; i++) {
     new_bins[i] = pow(10, from + i * width);

   }
   axis->Set(bins, new_bins);
   delete new_bins;
}


TH2D M_tau21_1("M_tau21_1", "M_{J} vs. #tau_{21} (#beta=1)",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau31_1("M_tau31_1", "M_{J} vs. #tau_{31} (#beta=1)",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau32_1("M_tau32_1", "M_{J} vs. #tau_{32} (#beta=1)",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau21_2("M_tau21_2", "M_{J} vs. #tau_{21} (#beta=2)",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau31_2("M_tau31_2", "M_{J} vs. #tau_{31} (#beta=2)",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau32_2("M_tau32_2", "M_{J} vs. #tau_{32} (#beta=2)",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau21_kt("M_tau21_kt", "M_{J} vs. #tau_{21} (kt)",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau31_kt("M_tau31_kt", "M_{J} vs. #tau_{31} (kt)",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau32_kt("M_tau32_kt", "M_{J} vs. #tau_{32} (kt)",
	     40, 0, 1, 20, 0, 1000);

TH2D M_tau21_prod("M_tau21_prod", "M_{J} vs. #Pi#tau_{21}",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau31_prod("M_tau31_prod", "M_{J} vs. #Pi#tau_{31}",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau32_prod("M_tau32_prod", "M_{J} vs. #Pi#tau_{32}",
	     40, 0, 1, 20, 0, 1000);

TH2D M_tau21_sum("M_tau21_sum", "M_{J} vs. #Sigma#tau_{2}/#Sigma#tau_{1}",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau31_sum("M_tau31_sum", "M_{J} vs. #Sigma#tau_{3}/#Sigma#tau_{1}",
	     40, 0, 1, 20, 0, 1000);
TH2D M_tau32_sum("M_tau32_sum", "M_{J} vs. #Sigma#tau_{3}/#Sigma#tau_{2}",
	     40, 0, 1, 20, 0, 1000);

TH1D M_cutoff("M_cutoff", "M_{J} cutoff weight", 50, 0, 1000);
TH1D M_original("M_original", "M_{J} original weight", 50, 0, 1000);
TH1D M_rescale("M_rescale", "M_{J} rescale weight", 50, 0, 1000);

TH1D M_original_unweight("M_original_unweight", "M_{J} original weight", 1, 0, 1000);
TH1D M_cutoff_unweight("M_rescale_unweight", "M_{J} rescale weight", 1, 0, 1000);


TH1D jet_pt1("jet_pt1", "Leading Jet p_{T}", 24, 0, 1200);
TH1D jet_pt2("jet_pt2", "Second Jet p_{T}", 20, 0, 1000);
TH1D jet_pt3("jet_pt3", "Third Jet p_{T}", 10, 0, 500);
TH1D jet_pt4("jet_pt4", "Fourth Jet p_{T}", 10, 0, 500);

TH1D jet_m1("jet_m1", "Lead Jet mass", 10, 0, 500);
TH1D jet_m2("jet_m2", "Second Jet mass", 10, 0, 500);
TH1D jet_m3("jet_m3", "Third Jet mass", 10, 0, 200);
TH1D jet_m4("jet_m4", "Fourth Jet mass", 10, 0, 200);



TH1D weight_bin1("weight_bin1", "Event weights for M_{J} #in (200, 400)", 100, -14, 6);
TH1D weight_bin2("weight_bin2", "Event weights for M_{J} #in (400, 600)", 100, -14, 6);
TH1D weight_bin3("weight_bin3", "Event weights for M_{J} #in (600, 800)", 100, -14, 6);
TH1D weight_bin4("weight_bin4", "Event weights for M_{J} #in (800, 1000)", 100, -14, 6);



double max_weight=1e99;
//double all_weight=0;
double reweight=1.0;
const double Mt_cutoff=400;
double weight_percentage=1;

void Ntuple::Begin(TTree * /*tree*/)
{
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).

  BinLogX(&weight_bin1);
  BinLogX(&weight_bin2);
  BinLogX(&weight_bin3);
  BinLogX(&weight_bin4);



  M_tau21_1.Sumw2();
  M_tau31_1.Sumw2();
  M_tau32_1.Sumw2();

  M_tau21_2.Sumw2();
  M_tau31_2.Sumw2();
  M_tau32_2.Sumw2();

  M_tau21_kt.Sumw2();
  M_tau31_kt.Sumw2();
  M_tau32_kt.Sumw2();

  M_tau21_prod.Sumw2();
  M_tau31_prod.Sumw2();
  M_tau32_prod.Sumw2();
  
  M_tau21_sum.Sumw2();
  M_tau31_sum.Sumw2();
  M_tau32_sum.Sumw2();

  jet_pt1.Sumw2();
  jet_pt2.Sumw2();
  jet_pt3.Sumw2();
  jet_pt4.Sumw2();
  
  jet_m1.Sumw2();
  jet_m2.Sumw2();
  jet_m3.Sumw2();
  jet_m4.Sumw2();

  TString option = GetOption();

}

void Ntuple::SlaveBegin(TTree * /*tree*/)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).
  TString option = GetOption();

  string opt(option);
  if(opt.find("_1.")!=string::npos)
    {weight_percentage=1e-2;}
  else if(opt.find("_2.")!=string::npos)
    {weight_percentage=2e-2;}
  else if(opt.find("_5.")!=string::npos)
    {weight_percentage=5e-2;}
  else if(opt.find("_8.")!=string::npos)
    {weight_percentage=8e-2;}
  else if(opt.find("_0.1.")!=string::npos)
    {weight_percentage=.1e-2;}
  else if(opt.find("_0.2.")!=string::npos)
    {weight_percentage=.2e-2;}
  else if(opt.find("_0.5.")!=string::npos)
    {weight_percentage=.5e-2;}
  else if(opt.find("_0.8.")!=string::npos)
    {weight_percentage=.8e-2;}
  else if(opt.find("_0.01.")!=string::npos)
    {weight_percentage=0.01e-2;}
  else if(opt.find("_0.02.")!=string::npos)
    {weight_percentage=0.02e-2;}
  else if(opt.find("_0.05.")!=string::npos)
    {weight_percentage=0.05e-2;}
  else if(opt.find("_0.08.")!=string::npos)
    {weight_percentage=0.08e-2;}
  else if(opt.find("_0.001.")!=string::npos)
    {weight_percentage=.001e-2;}
  else if(opt.find("_0.002.")!=string::npos)
    {weight_percentage=.002e-2;}
  else if(opt.find("_0.005.")!=string::npos)
    {weight_percentage=.005e-2;}
  else if(opt.find("_0.008.")!=string::npos)
    {weight_percentage=.008e-2;}
  else if(opt.find("_0.0001.")!=string::npos)
    {weight_percentage=.0001e-2;}
  else if(opt.find("_0.0002.")!=string::npos)
    {weight_percentage=.0002e-2;}
  else if(opt.find("_0.0005.")!=string::npos)
    {weight_percentage=.0005e-2;}
  else if(opt.find("_0.0008.")!=string::npos)
    {weight_percentage=.0008e-2;}


  cout<<"Message: Weights cutoff "<<weight_percentage*100<<"%"<<endl;


}

Bool_t Ntuple::Process(Long64_t entry)
{
   // The Process() function is called for each entry in the tree (or possibly
   // keyed object in the case of PROOF) to be processed. The entry argument
   // specifies which entry in the currently loaded tree is to be processed.
   // It can be passed to either Ntuple::GetEntry() or TBranch::GetEntry()
   // to read either all or the required parts of the data. When processing
   // keyed objects with PROOF, the object is already loaded and is available
   // via the fObject pointer.
   //
   // This function should contain the "body" of the analysis. It can contain
   // simple or elaborate selection criteria, run algorithms on the data
   // of the event and typically fill histograms.
   //
   // The processing can be stopped by calling Abort().
   //
   // Use fStatus to set the return value of TTree::Process().
   //
   // The return value is currently not used.

  //obtain M_J dependent scaling factor


  if(entry==0 && b_QCD)
    {
      max_weight=0;
      
      for(Long64_t i=0; i<fChain->GetEntries(); i++)
	{
	  fChain->GetEntry(i);
	  weight*=scale;
	  
	  double Mt=0;
	  for(int j=0; j<n ;j++)
	    if(pt[j]>50.0)
	      Mt+=m[j];
	  
	  if(Mt>400)
	    max_weight+= weight;
	}  
      /*
      double cut_weight=0;
      for(Long64_t i=0; i<fChain->GetEntries(); i++)
	{
	  fChain->GetEntry(i);
	  weight*=scale;
	  if(weight > weight_percentage*max_weight)
	    continue;
	  cut_weight+= weight;
	}
      */
      reweight=1;//max_weight/cut_weight;
    }

  fChain->GetEntry(entry);
  if(n==0) return kTRUE;

  if(!b_QCD)
    weight=1.0;

  //scale according to cross sections
  weight*=scale;
 
  double Mt=0;
  for(int i=0; i<n ;i++)
    if(pt[i]>50.0)
      Mt+=m[i];
  
  if(Mt>400 && b_QCD)
    {
      if(weight > weight_percentage*max_weight)
	return kTRUE;
    }
  
  double tau1_sum=0;
  double tau2_sum=0;
  double tau3_sum=0;
  double tau1_prod=1;
  double tau2_prod=1;
  double tau3_prod=1;
  

  for(int i=0; i<n ;i++)
    if(pt[i]>50.0)
      {
	if(i<4)
	  {
	    tau1_sum+=tau1_1[i];
	    tau2_sum+=tau2_1[i];
	    tau3_sum+=tau3_1[i];

	    tau1_prod*=tau1_1[i];
	    tau2_prod*=tau2_1[i];
	    tau3_prod*=tau3_1[i];
	  }
      }

  if(Mt>=200)
    {
      if(Mt<400)
	weight_bin1.Fill(weight);
      else if(Mt<600)
	weight_bin2.Fill(weight);
      else if(Mt<800)
	weight_bin3.Fill(weight);
      else
	weight_bin4.Fill(weight);
    }

  if(Mt>=400)
    {
      jet_pt1.Fill(pt[0], weight);
      jet_pt2.Fill(pt[1], weight);
      jet_pt3.Fill(pt[2], weight);
      jet_pt4.Fill(pt[3], weight);
      
      jet_m1.Fill(m[0], weight);
      jet_m2.Fill(m[1], weight);
      jet_m3.Fill(m[2], weight);
      jet_m4.Fill(m[3], weight);

    }
    


  /*
  if((b_QCD && Mt>=Mt_cutoff && weight >= max_weight))
    return kTRUE;
  */
   
  double weight_Mt=weight;

  /*
  int c_bin=M_rescale.FindBin(Mt);  
  if(b_QCD && Mt>=Mt_cutoff)
    weight_Mt=weight*M_rescale.GetBinContent(c_bin);
  else
    weight_Mt=weight;
  */

  M_tau21_1.Fill(tau2_1[0]/tau1_1[0], Mt, weight_Mt);
  M_tau31_1.Fill(tau3_1[0]/tau1_1[0], Mt, weight_Mt);
  M_tau32_1.Fill(tau3_1[0]/tau2_1[0], Mt, weight_Mt);
  
  M_tau21_2.Fill(tau2_2[0]/tau1_2[0], Mt, weight_Mt);
  M_tau31_2.Fill(tau3_2[0]/tau1_2[0], Mt, weight_Mt);
  M_tau32_2.Fill(tau3_2[0]/tau2_2[0], Mt, weight_Mt);

  M_tau21_kt.Fill(tau2_kt[0]/tau1_kt[0], Mt, weight_Mt);
  M_tau31_kt.Fill(tau3_kt[0]/tau1_kt[0], Mt, weight_Mt);
  M_tau32_kt.Fill(tau3_kt[0]/tau2_kt[0], Mt, weight_Mt);

  M_tau21_prod.Fill(pow(tau2_prod/tau1_prod,0.25), Mt, weight_Mt);
  M_tau31_prod.Fill(pow(tau3_prod/tau1_prod,0.25), Mt, weight_Mt);
  M_tau32_prod.Fill(pow(tau3_prod/tau2_prod,0.25), Mt, weight_Mt);
  
  M_tau21_sum.Fill(tau2_sum/tau1_sum, Mt, weight_Mt);
  M_tau31_sum.Fill(tau3_sum/tau1_sum, Mt, weight_Mt);
  M_tau32_sum.Fill(tau3_sum/tau2_sum, Mt, weight_Mt);


  return kTRUE;
}

void Ntuple::SlaveTerminate()
{
   // The SlaveTerminate() function is called after all entries or objects
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server.

}

void Ntuple::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
  TFile fileout(GetOption(), "RECREATE");


  M_tau21_1.SetOption("colz");
  M_tau31_1.SetOption("colz");
  M_tau32_1.SetOption("colz");

  M_tau21_2.SetOption("colz");
  M_tau31_2.SetOption("colz");
  M_tau32_2.SetOption("colz");

  M_tau21_kt.SetOption("colz");
  M_tau31_kt.SetOption("colz");
  M_tau32_kt.SetOption("colz");

  M_tau21_prod.SetOption("colz");
  M_tau31_prod.SetOption("colz");
  M_tau32_prod.SetOption("colz");
  
  M_tau21_sum.SetOption("colz");
  M_tau31_sum.SetOption("colz");
  M_tau32_sum.SetOption("colz");



  
  M_tau21_1.Write();
  M_tau31_1.Write();
  M_tau32_1.Write();

  M_tau21_2.Write();
  M_tau31_2.Write();
  M_tau32_2.Write();

  M_tau21_kt.Write();
  M_tau31_kt.Write();
  M_tau32_kt.Write();

  M_tau21_prod.Write();
  M_tau31_prod.Write();
  M_tau32_prod.Write();
  
  M_tau21_sum.Write();
  M_tau31_sum.Write();
  M_tau32_sum.Write();


  M_original.Write();
  M_cutoff.Write();
  M_original_unweight.Write();
  M_cutoff_unweight.Write();
  M_rescale.Write();

  weight_bin1.Write();
  weight_bin2.Write();
  weight_bin3.Write();
  weight_bin4.Write();

  jet_pt1.Write();
  jet_pt2.Write();
  jet_pt3.Write();
  jet_pt4.Write();
  
  jet_m1.Write();
  jet_m2.Write();
  jet_m3.Write();
  jet_m4.Write();
  
  
  fileout.Close();

}
